home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / rexx / imc / rexx-imc.5 / rxmathfn.exec < prev    next >
Encoding:
Text File  |  1993-06-25  |  6.6 KB  |  214 lines
  1. /* Example Rexx program to provide mathematical functions */
  2. trace off
  3.  
  4. Ecall=40 /* Incorrect call to routine */
  5. Enum=41  /* Bad arithmetic conversion */
  6. Eundef=43   /* Routine not found */
  7.  
  8. parse upper source . . . what ./* Must have what equalling the function name */
  9. if Arg()<>1+(what="TOPOWER") then do
  10.    say what':'errortext(Ecall)
  11.    return   /* If an error occurs, return with no value causes the caller */
  12. end         /* to be flagged (with "Function did not return result")      */
  13.  
  14. parse arg in .,in2 .
  15. if \datatype(in,'n') then do
  16.    say what':'errortext(Enum)
  17.    return
  18. end
  19.  
  20. d=digits()
  21. f=fuzz()
  22.  
  23. if d<20 then d1=d;else d1=19
  24. lim=10**-d1/20
  25. numeric digits d1+2
  26. pi=3.1415926535897932384626433
  27. halfpi=pi/2
  28. twopi=pi+pi
  29. log10=2.3025850929940456840179910
  30. log2=0.6931471805599453094172321
  31. log10two=0.30102999566398119521373889
  32.  
  33. if d<20 then d1=d;else d1=19
  34. lim=(1'E-'d1)/20
  35.  
  36. select
  37.    when what="SIN" then answer = sin(in)
  38.    when what="COS" then answer = sin(in+halfpi)
  39.    when what="TAN" then answer = sin(in)/sin(in+halfpi)
  40.    when what="EXP" then answer = exp(in)
  41.    when what="ATAN" then answer = atn(in)
  42.    when what="LN" then answer = ln(in)
  43.    when what="TOPOWER" then do
  44.          if in2=0 then answer = 1
  45.          else if in=0 & in2>0 then answer = 0
  46.          else answer = exp(ln(in)*in2)
  47.    end
  48.    when what="ASIN" then answer = 2*atn(in/(1+sqr1(1-in*in)))
  49.    when what="ACOS" then answer = halfpi-2*atn(in/(1+sqr1(1-in*in)))
  50.    when what="SQRT" then return sqr(in)
  51.    otherwise say what':'errortext(Eundef) ; return
  52. end
  53. numeric digits d1
  54. answer=format(answer)
  55. return answer
  56.  
  57. series:     /* Calculates a Chebyshev series with the given coefficients, */
  58. n=arg()-1   /* evaluated at the given argument.  arg(1) is the argument, and */
  59. parse arg z /* all the other arg() are the series coefficients */
  60. z=z+z;m2=0;t=0
  61. do i=n+1 to 2 by -1
  62.    m1=m2
  63.    u=t*z-m2+arg(i)
  64.    m2=t
  65.    t=u
  66. end
  67. return t-m1
  68.  
  69. sin: parse arg x
  70. /* reduce range of x */
  71. x=(x/halfpi+2)//4-2
  72. if x<-1 then x=-x-2
  73. if x>1 then x=2-x
  74.  
  75. return series(2*x*x-1,1.27627896240226588021,-0.14263078459051800478,,
  76.    0.00455900800332590124,-0.00006829375677098333,0.00000059248092883084,,
  77.    -0.00000000335139580191,0.00000000001333639299,-0.00000000000003936461,,
  78.    0.00000000000000008961,-0.00000000000000000016),
  79.  *x
  80.  
  81. exp: parse arg x
  82. x=x/log10
  83. e=trunc(x)
  84. x=x-e
  85. if x<0 then do;x=x+1;e=e-1;end
  86. /* now have answer = 10**x * 10**e where 0<=x<1 and e is an integer */
  87. return series(x+x-1,4.30022915484994695129,2.13908204064962629380,,
  88.      0.58426304848002648400,0.10914429717084543695,0.01545385656994400473,,
  89.      0.00175990305162950212,0.00016753288979627516,0.00001369642180963234,,
  90.      0.00000098103821148018,0.00000006251839260577,0.00000000358806613888,,
  91.      0.00000000018729961298,0.00000000000896585002,0.00000000000039629176,,
  92.      0.00000000000001626892,0.00000000000000062348,0.00000000000000002240,,
  93.      0.00000000000000000075,0.00000000000000000002)'E'e
  94.  
  95.  
  96. atn: parse arg x
  97. select
  98.    when x<-1 then do;w=-halfpi;y=-1/x;end
  99.    when x>1 then do;w=halfpi;y=-1/x;end
  100.    otherwise w=0;y=x
  101. end
  102. /* now -1<=y<=1 and atn x=w+atn y */
  103. return series(2*y*y-1,0.88137358701954302523,-0.05294646227335292762,,
  104.    0.00556792102970276495,-0.00069059750180019885,0.00009287148663927100,,
  105.    -0.00001310759805636663,0.00000191051829724472,-0.00000028495930832883,,
  106.    0.00000004324438932225,-0.00000000665169198781,0.00000000103425288179,,
  107.    -0.00000000016224319604,0.00000000002563983127,-0.00000000000407739378,,
  108.    0.00000000000065190359,-0.00000000000010471489,0.00000000000001688910,,
  109.    -0.00000000000000273382,0.00000000000000044394,-0.00000000000000007229,,
  110.    0.00000000000000001180,-0.00000000000000000193,0.00000000000000000031,,
  111.    -0.00000000000000000005)*y+w
  112.  
  113. ln: parse arg x
  114. if x<=0 then do
  115.    say 'ln('||x'):'errortext(Ecall)
  116.    return
  117. end
  118. if x=1 then return 0
  119. if x<1 then do;s='-';x=1/x;end
  120. else s=''
  121. numeric form scientific
  122. parse value format(x,1,,,0) with . 'E' exp
  123. if exp='' then exp=0
  124. e=trunc((1+exp)/log10two-3)
  125. /* x/(2**e) is now between 1 and 16 */
  126. x=x/2**e
  127. do while x>1.6
  128.    x=x/2
  129.    e=e+1
  130. end
  131. ans=series(2.5*x-3,0.93022922133637741858,-0.08184145665572451963,
  132.    ,0.00947661158839661302,-0.00122828373982062934,0.00016939529033498945,
  133.    ,-0.00002430127203947124,0.00000358276167843137,-0.00000053890849196729,
  134.    ,0.00000008231521903733,-0.00000001272671619709,0.00000000198712378075,
  135.    ,-0.00000000031280003537,0.00000000004957693767,-0.00000000000790361608,
  136.    ,0.00000000000126636133,-0.00000000000020379548,0.00000000000003292356,
  137.    ,-0.00000000000000533708,0.00000000000000086781,-0.00000000000000014149,
  138.    ,0.00000000000000002312,-0.00000000000000000378,0.00000000000000000062,
  139.    ,-0.00000000000000000010,0.00000000000000000001)*(x-1)+e*log2
  140. return s||ans
  141.  
  142. sqr1: /* just return exp(ln(x)/2) unless x=0 */
  143. parse arg z
  144. if z=0 then return 0
  145. return exp(ln(z)/2)
  146.  
  147. sqr2: /* use digit-by-digit algorithm to find sqr(x) to required number of
  148.          places.  This routine is no longer in use. */
  149. parse arg a
  150. numeric digits d+1
  151. if a=0 then b=0
  152. else if a<0 then do
  153.    say "sqrt("a"):"errortext(Ecall)
  154.    return
  155. end
  156. else do
  157.    numeric form scientific
  158.    parse value format(a,1,,,0) with before '.' after 'E' e
  159.    if pos('E',before)>0 then parse var before before 'E' e
  160.    a=before||after  /* a is now pointless */
  161.    if e='' then e=0 /* e is the exponent  */
  162.    if e//2 \= 0 then e=e-1
  163.    else a='0'a      /* exponent is odd, and decimal comes after 2 digits */
  164.    b='';b1=0;c=0
  165.    numeric digits d%2+2 /* avoid unnecessary work - we only need this precision*/
  166.    do i=1 by 2 until length(b)=d+1
  167.       b1=trunc(b1)||substr(a,i,2,0)
  168.       /* find least x s.t. (c||x)*x>b1. Find (c||x-1)*(x-1). */
  169.       d2=c||1
  170.       d1=0
  171.       x=1
  172.       do while d2<=b1
  173.          d1=d2
  174.          d2=d2+(c||x)
  175.          x=x+1
  176.          d2=d2+x
  177.       end
  178.       x=x-1
  179.       b=b||x
  180.       b1=b1-d1
  181.       c=trunc((c||x)+x) /* never in exponential form */
  182.    end
  183.    b=left(b,1)'.'substr(b,2)'E'e/2
  184. end
  185. numeric digits d
  186. return format(b)
  187.  
  188. sqr: /* Use Newton's method to get square root to required number of places */
  189. parse arg x
  190. if x=0 then return 0
  191. if x<0 then do
  192.    say "sqrt("||x"):"errortext(Ecall)
  193.    return
  194. end
  195. numeric form scientific
  196. parse value format(x,1,1,,0) with n +1 'E' exp /* error if x<0 */
  197. if exp='' then exp=0
  198. r=n/2||'E'||(exp+1)%2 /* divide both mantissa and exponent by 2 for initial
  199.                          estimate */
  200. d1=5
  201. numeric digits d1+2
  202. numeric fuzz 1
  203. do until r = r1       /* Get estimate accurate up to d1 digits */
  204.    parse value r (x/r+r)/2  with r1 r
  205. end
  206. do while d1<d
  207.    d1=min(d1+d1,d)    /* Double up until required accuracy is reached */
  208.    numeric digits d1+2
  209.    parse value r (x/r+r)/2  with r1 r
  210. end
  211. numeric digits d
  212. return format(r)
  213.  
  214.